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I. INTRODUCTION 


This report presents the results of a study on the evaluation and 
reduction of noise impact to a community due to aircraft operation. 
Existing techniques have been used to assess the noise impact and to 
optimize the flight paths of an approaching aircraft with respect to the 
annoyance produced. Major achievements have been! (1) the development 
of a population model suitable for determining the noise impact, 

(2) generation of a numerical computer code which uses this population 
model along with the steepest descent algorithm to optimize approach/ 
landing trajectories, (3) implementation of this optimization code 
in several fictitious cases as well as for the community surrounding 
Patrick Henry International Airport. 

Previous work has centered on developing noise annoyance criteria 
for flyover (i.e. NEF, NNI, CNR, etc.) and ground noise signatures for 
aircraft. Some of these criteria are discussed in References 1-5 with 
a review of many of the noise effect measures being summarized in Ref. 6. 
Typical of the noise footprint work is Ref. 7. The annoyance criterion 
used in the study is the noise impact index (Nil) . 

The details of the models used, their advantages and disadvantages 
and the results obtained are outlined in the following sections. 
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II. PROBLEM FORMULATION 


A. Overview 

Analysis of the problem consists of six parts: (1) aircraft 

noise signatures, (2) population models, (3) cost (annoyance) function, 
(4) aircraft flight-path model, (5) aircraft constraints, and 
(6) approach/landing path optimization. A modular concept has been 
employed so that modification of any of these segments may be effected 
with relative ease. The sections below describe each of these parts 
in detail. 

B. A/C Noise Signature 

The aircraft noise signature is obtained using data from Ref. 3. 
Here the effective perceived noise level (EPNdB) is given as a function 
of slant range to the closest point of approach for a variety of 
aircraft. A typical plot of the slant range variation is shown in 
Figure 1. These data were fit using standard least squares techniques 
to yield an expression for EPNdB given by 

EPNdB = 115-22.5 log^o x ( slant Range). (D 

This equation is used for calculation of the maximum noise level at 
each location for a flyover. A typical footprint for a straight in 
approach along a 3-degree glide slope is shown in Figure 2. 

C. Population Model - 

To model the population, a map of the community is overlaid with 

a grid and the population in each section of the grid determined. The 
population distribution within each section is assumed to be uniform. 
Several grid geometries were examined (see Figure 3). These geometries 
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3. Concentric Circles 


Figure 3. Grid Geometries. 
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included: (1) rectangular sections of equal size, (2) rectangular 

sections whose dimensions increased with distance from the airport 
runway, and (3) concentric circles divided by several radial lines. 

The second scheme was chosen since it requires fewer rectangular 
sections than the first and is easier to implement than the third. 
Computer time required for determining the optimum trajectory varies 
directly with the number of grid sections. Furthermore, in light of 
the dependence of noise levels on distance and the fact that the 
aircraft has higher altitude when further from the runway, the need for 
high resolution of the population density diminishes with distance from 
the airport. While the third population scheme could offer this same 
advantage, it is somewhat more difficult to determine the population 
and noise impact for each grid section with such a geometry. 

Within a grid section, the population is determined by use of the 
SITE II system, (Ref. 8), available on the CDC 7600 computer at the 
NASA-Langley facility. This system requires as input the latitude and 
longitude of a reference point and the coordinates of the corners of 
each rectangular section. Although SITE II allows for simple retrieval 
of 1970 census data, there is some question about its resolution 
capabilities for small grid sections. In addition, in rapidly growing 
areas the population data may lag actual population. The SITE II 
program is capable of producing detailed census information as shown 
in Figure A. However for the present analysis only population 
information is used. 
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Figure 4. Demographic Profile Report 
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D. Flight Path Model 


There are two ways in which the trajectory of the aircraft can be 
determined. In one, a discrete time integration of the equations of 
motion with control deflections yields point by point spatial coordinates 
and orientation. Although this allows the flexibility of building in 
control constraints as well as dynamical constraints (e.g. max roll angle) 
it requires a considerable number of states to be stored in the 
optimization routine. In a multi-aircraft, multi-runway problem, as is 
anticipated, these storage requirements become prohibitive. 

Thus, another method was adopted which utilizes only the functional 
form of the trajectory to describe the flight path. First a starting 


path was assumed which went from the initial point to the desired 
runway and ended up with the proper heading, i.e., velocity vector 
aligned with the runway. The following equation was used to generate 
this starting trajectory. (See Figures). 


y s ( x ) 



(x-Xp) + (y p -y 0 ) 


EXP j-C(x-Xf) / (x Q -xf )J + y 0 

( 2 ) 


For the vertical motion a simple three degree descent path was 
assumed. 

Next the first five Fourier sine harmonics were used to introduce 


deviation from this starting path. One advantage in using this type 
representation is the fact that each of the forms contributes zero at 
the end points. Therefore if the starting path satisfies the boundary 
conditions the curve with the deviations will also. An exponential 
decay at the final point was used to eliminate heading deviations. 
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The equations with the deviations thus become 


y(x) 


a i sin[iTi(x-x 0 )/(x f -x 0 )]^| |l-EXP[(x-Xf)/Ci]^| + y s (x) 

(3a) 


Z(x) =< £ BiSin[iri (x-x Q ) / (x f -x a ) ]\ J l-EXP[(x-Xf)/Ci] 1+ Z g (x) 

V 1 A J (3b) 

A second advantage to using a Fourier Series representation for the 


curve is that it provides a means of representing a function with a finite 
number of parameters. This reduces the optimization problem from a 
variational one to an ordinary one. 

E. Constraints 

The use of a functional form of the flight path for the trajectory 
requires the reformulation of constraints into parameters which can be 
used in the optimization. This is accomplished by translating the 
steady state solutions of the lateral and longitudinal perturbation 
equations into geometric constraints. An exact derivation is given in 
Appendix A. The constraints are incorporated by determining maximum 
curvature and slope parameters as a function of aerodynamics and 
physical constraints. For example the constraint of a maximum roll 
angle, <J> max , yields 


d 2 y 

— — *T C 0 +C.Co (J) 

dx^ < 2 Id 't'max _ ^ 

1+ / dy_\ 2 - C 4 +C 1 C 5 V avg 

Mx' 

where Cq through C5 depend upon aircraft stability and control derivations 
(see Appendix A for details) and V avg is the average velocity. Similar 
expressions are given in Appendix A for constraints on aileron rudder 
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and elevator deflection, flight path angle and pitch rate limits. 

F. Cost Function 

A large number of criteria have been proposed by evaluating noise 
annoyance (e.g., EPNdB, Nil, sleep interference index, speech inter- 
ference index, etc.). The recent trend in noise assessment work is 
toward a universal measure — the noise impact index (Nil) . This measure 
is a weighted day— night model which accounts for population density. 

It is described in detail in Ref. 9. Briefly, the total population 
exposed to each incremental average day-night model sound level is 
multiplied by the weighting function for the level. The weighting 
function used in shown in Figure 6. This weighting factor W(L dn ) 
multiplied by the population exposed to that L dn is summed and 
normalized by the total population giving the Noise Impact Index for 


l P(L dn )W(L dn ) 
L dn 

I p a dn ) 

^dn 


(5) 


The cost function or payoff for the optimization procedure is 
taken to be the Nil plus penalties for violating constraints. Basically 
the optimization procedure is set up to "drive” the aircraft trajectory 
to the path which will minimize the Nil and at the same time not 
violate any constraints. As an example of the constraint of flight 
path angle not exceeding a maximum descent angle, Y d , nor a maximum 
climb angle, Y c , is written as 

tanY c < < tanY d (6) 
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SOUND LEVEL WEIGHTING FUNCTION 


FOR OVERALL IMPACT ANALYSIS 



Figure 6. Sound Level Weighting Function for Overall 
Impact Analysis. 


12 




Each is converted to a penalty which is added to the Nil in the form 

Cost = NII+(||/tanY d ) 2Q + (tanY c /||) 20 (7) 

As is seen for values of the flight path angle within the allowable 
range the penalty is negligible; however for values outside this range 
the penalty and thus the increase in cost is great. Other terms are 
added in a like manner. 
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III. OPTIMIZATION 


The opt imum trajectory is determined by calculating values of the 
cii's and B-j/s (Eq. 2) which minimize the total cost (Nil plus penalties). 
A steepest descent algorithm is employed here. Basically, this method 
computes the gradient of the cost function, C, with respect to the s 
and Bi's, then searches along the negative gradient direction for values 


of a^s and Bi's which reduce the cost. 


5 

AC = l 
i=l 


,3C . , 3C A0 n 

( 3ST + 


The change in cost is given by 

( 8 ) 


The process continues iteratively until the cost converges to within a 
specified tolerance. While implementation of the algorithm is fairly 
straightforward, convergence near the optimal set of cx^'s and s is 
inherently slow. Most of the cost reduction, however, occurs in the 
first few iterations. 

A. The Optimization Algorithm 

A computer code has been developed which implements the functions 
described above. Figure 7 shows a flow chart for this code. Initial 
data (population map, aircraft constraints, initial and final aircraft 
positions, etc.) are required for each airport/airplane configuration 
to be evaluated. To facilitate calculation of the Fourier coefficients, 
the coordinate axes are rotated such that a line joining the initial and 
final aircraft positions is made to be parallel to the x axis. A nominal 
trajectory is generated which constrains the heading of the aircraft to 
asymptotically approach the runway. The steepest descent search then 
begins and continues until the stopping criterion is met. This criterion 
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is met if successive improvements become negligible. 

In order to provide a more accurate noise impact in each population 
section the impact is integrated using quadratures. This procedure can 
be found in Ref. 9. 

The various functions such as the population model, the cost function, 
and the aircraft signature are incorporated as subroutines. This will 
allow ease of upgrading or modification if different models are desired. 

Appendix B contains the Fortran code as written for a CDC Cyber 172 
machine. 

B. Results 

Several cases have been run to test the benefits that can be 
obtained by this approach. First, a fictitious set of data incorporating 
a population valley is used. As can be seen in Figure 8 the optimization 
algorithm moves the aircraft (a Convair 880) towards the valley (i.e. 
fewer people impacted) with a corresponding improvement in the Nil of 
32%. 

The second case models the Patrick Henry Airport in Hampton, 

Virginia. Here the SITE II program was used to generate the census 
data for each block as shown in Figure 9. Two initial trajectories 
were flown. One entering the area from the northwest over the 
Swing VOR station and the other from the southwest over the Franklin 
VOR station. Both of these paths are specified IFR trajectories 
(Figure 10). The aircraft enters the area approximately 30,000 meters 
from the runway. In addition several straight-in paths were 
evaluated. Figure 9 shows each of the trajectories. The associated 
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Figure 3. Optimization Results Using Fictitious Population Data. 



Figure 9. Population Model and Optimization Results for Patrick Henry Airport. 
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Table I 


Northwest Approach 
Entry Point: Swing 


Traj. 

No. 

Description 

Cost 

(Nil x 10 -2 ) 

% Change from 
Present 





1 

60 deg wrt runway 

2.373 

+3.2% 

2 

30 deg wrt runway 

2.438 

+6.0% 

3 

| 

Initial iteration 

2.27 

-1.3% 

! 

4 

Optimal 

2.213 

-3.8% 

5 

Straight in 

2.316 

+.1% 

N1 

Presently used 

2.300 

0% 


Table II 

Southwest Approach 
Entry Point : Franklin 


Traj . 
No. 

Description 

Cost 

(Nil x 10 -2 ) 

% Change From 
Present 

5 

Straight in 

2.316 

-1.3% 

6 

Initial iteration 

2.408 

+2.6% 

7 

Optimal 

2.241 

- 4.5% 

8 

30 deg wrt runway 

2.598 

+10.7% 


60 deg wrt runway 

2.687 

+14.5% 


Presently used 

2.346 

0% 
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Nil's are summarized in Tables I and II. 

As is seen that even for this case, where the population is sparse 
away from the runway and congested near the end of it, an improvement 
of 3 to 5% is achieved using the optimization algorithm. It should 
also be noted that this technique allows not only the optimization of 
the path but also can be used to evaluate existing or proposed paths, 
such as the nominal and straight-in paths indicated. 

Conclusion 

A method has been formulated to optimize the path of an aircraft 
during approach or take-off from any airport. Models have been developed 
using available data where possible for population, aircraft signature, 
noise impact, constraints and flight path. An algorithm using steepest 
descent has been implemented and tested. This approach allows 

1) The evaluation of the noise impact of existing flight paths, 

2) The evaluation of the noise impact of proposed flight paths, and 

3) The optimization of the flight path to minimize the noise 
impact under constraints. 

This method has been applied to the Patrick Henry International Airport. 
Both nominal and other straight paths were evaluated. Also an optimal 
path was determined for each of two terminal area entry points. 
Performance ranged from a 15% degredation of the Nil to 4.5% improvement 
compared to the presently used approaches. It is significant that as 
much as 3 to 5% improvement could be achieved in light of the fact that 
most of the population is concentrated at the end of the runway. 
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APPENDIX A 


Derivation of Parameterized Trajectory Constraints 
Lateral perturbation equations 
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Taking the Laplace transform (I.C.’s - 0) 

L eq'n: (i x s 2 -C. s) + (s) + (-i x Zs 2 -C s) Y (a) = C. 6 a (s)+C^ « r (s) 
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A-l 
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To determine the required <5 a for a given 6 r we consider 6 a an unknown 
along with <J>(s) and 'P(s) [i.e. move to the left hand side of the 
equations] and solve for <5 a /6 r using Cramer's rule 
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The denominator (characteristic eqn.) is given by: 

Ms) = s 4 { -Cy 6 a (i x i z -i2z) }+s 3 { c y 5 a [ i zV i Ar+ ix z C^n p ) ] 
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The numerator is: 
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( 6 ) 


Now assuming that only the steady state (st. st.) condition is of 

interest. 

Lira N(s) = ,6^. 

s+0 A (s) ^6/ St ' St * 


we get 


6 a 

(7—) st. st. = 
°r 


~ c ng r (g2 G £p~gi c 2. r )~ c g,iS r ^g2 Cn p g l Cn r^ 

Cn 6 a )+Cs '6a ^ g 2 TJ n p " g i T: n r ^ 

cos0 o (C n6 . C^+C^ C nr )-sin0 o (C n ^ C £p + ^’£ ( $ ^np) 


(7) 


(^) st. st. = -cos0 o (C n6a C £r +C A6a C nr )+sin0 o (C n ^C )lp +C z ^C np ) (8) 
For small initial flight path angle (i.e. 0 o ~O) 

c n 6r c £ r +C £5 r c n r 


(57) st * st ‘ = " SrV 1 *^ " C 1 _ 


: C n r “ C ] (9) 

Assuming 0 o =O to simplify we can write the transfer functions for <|> and 







— 


c £ 5 C ne _C n5 % 



Y 

5 

— = ~ — = c 2 

(10) 

— 

r 

C £3 C n r -Cn e C Jlr 


— 

¥ 

<5„ 

% Cn e" c n 6 % 

■ 6 ; l a =c 3 

(ID 

— 

u a 

c £g c n r c ngC£ r 



4> _ 

C y$ ^ C ^ r C ng _c ^ 3 C n r ^ +c )l (Sr ^ c y p C n r +C n 3 ( m_c y r ) ) +c n 5 ^ ( c £g 

(m-Cy r )+CygC^ r ) 


6 r " 

^ C ^$ Cn r~ C ng C ^r ) 

n oo 



= 

c 4 

(12) 

— 

4 ) 

C y g ^ G Jl r c ng -C 2 ,g Cn r^ +C 2 . 6 a ^ C yg C nr +Cn g ^ m-C yr^ +c n 5a ( c s,g 

(m-Cyr ) +Cy g C ) 


The slope at any point is and the angle the slope makes with the 
x axis is tan 1 (^) . 

The angular rate 7 is then ^ tan -1 (^) 


or 


a_ 

9x 


{tan-1 & = V — { 

itan dx' dt V avg 3 x 




d 2 y 

f " (x) 

Then T - V avg _dxL - V avg < 1+ FTOOF’ (14) 


If we know 5 r we can determine S a from 6 a = Cj _6 r 
Also Y = C 2 6 r +C 3 5 a = (C 2 +C 1 +C 3 )5 r 


We can also write 

<J> = C45 r +C56 a = (C4+C2C5)6 r 


Constraining <5 a to be < 5 a max 

5- to be < 5 r 
r r max 

and <j> to be £ <J> 


max 


we get the following expressions 

<|>max 

r l - C 4 +C 2 C 5 


5 t* o 5 -w* 
r 2 r max 


6 r 3 < a max 




(~ max bank angle) 


(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 
( 21 ) 


A- 5 


( 22 ) 



The constraining value is given by 


6 r max " min ( 5 ri> 6 r 2 ’ 6 r 3 ) 


(23) 


which yields 

\iax = ( C 2 +C 1 C 3^ niin(6 r;L ,S r2 ,6 r3 ) 

This condition incorporates all three constraints ((17)-(19)) as 

d 2 

t 

1+ f'(x) 2 


(24) 


d^y 




e . mln (5ri>6r2 , Sr3 ) 


avg 


Longitudinally we wish to constrain the behavior of the trajectory so 
that we restrict y (the flight path angle) and 0 (the pitching rate) . 
The trajectory is given by 



Then, assuming the aircraft center of mass follows this trajectory Y is 
given by 


y 


tan' 


-1 


dz 

dx 



We wish to constrain y to a maximum descent angle, Yd ma x 3nd 3 maximuin 

an § le ’ Y <W 


dz 


tanY c max < s - tan Yd ®« 


A- 6 


Thus 
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5 


10 


15 


20 


25 


td 

H 30 


35 


40 


45 


50 


55 


PROGRAM NOISE 1 1 NPUT * OUTPUT . T A PE 5 = I NPUT t T APE6=0UTPUT * TAPE97 * TAPE9B 
1*TAPE99> 

COMMON ALF A ( 5 ) * BET A ( 5 ) * POS I T < 53 . 3 ) . ARR AY ( 578 , 9 ) ♦ NM AP 
COMMON /CURVE/ YCURVE ( 5 1 ) « ADT ( 51 ) , ADDY ( 51 ) 

COMMON /LABEL/ L INfO ( 4 J * LLOC ( 3 ) 

COMMON /AIRPORT/ XPORT * YPOR T * ZPORT 
COMMON /SCALE/ XM I N * X I NC » YMI N , Y I NC 


INTEGER CO JNT * HALF 

DIMENSION ALFA0D(5)* BETAOD ( 5) * GY(5)« GZ(5)« 0ALFM5)* 
DIMENSION AG Y ( 5 ) . 0GY(5)* AGZ(5)* RGZ(5) 


DBETA ( 5) 


C 

C 

C . 

C , READ MAP FROM DISC INTO MEMORY 
C . 

c 

c 

READ (5**) AU.A12 

READ (5**) NMAPiXPORTOi YPORTO 

READ ( 5 * * ) ( (ARRAY f I ♦ J) * J = 1 * 9 ) « I=1*NMAP) 

C 

C 

C . 

C . INPUT INITIAL CONDITIONS # 

C . 

C 

C READ (5.*) MAXIT*YALLOW*ZALLOW* ( ALFA( I ) *1 = 1*5) * (BETA( I) *1 = 1*5) * XU. 

IYOiZO.XF. YF.ZF 

READ (5*910 0) (LLOC(I) . 1 = 1*3) * (HNFO( I ) • 1 = 1*4) 

WRITE (6*9110) (LLnC(I) *1=1*3) * (LINFO(I) *1=1*4) 

WRITE (6.9120) MAXIT 

WRITE (6*9010) XO*YO*ZO*XF* YF * ZF » XPORT 0 » YpOKTO 
WRITE (6*9130) Y ALLOW * Z ALLOW 
WRITE (6*9020) A11.A12 

WRITE (6.9140) ( I ♦ ALF A ( I ) , BETA (I ) . 1 = 1 * 5 ) 


C 

C 

C 

C 

C 

C 

c 


! COORDINATE ROTATION 

t . 


THETA = ATAN2((YF-Y0)*(XF-X0)) 

A = ( XF*COS ( THETA ) + YF*SIfJ (THETA ))•( XO* COS ( THETA ) + Y0 + SIN (THETA) ) 

PHI = ATAN2( (ZF-ZO; .A) 

XOC AP = XO*COS(THETA)*COS(PHl ) ♦YO*COS ( PHI )*SIN( THETA )+ZO*SIN( PHI ) 
XFCAP = XF+COS ( THETA ) *C OS (PHI )+YF*SlN( THETA )*COS( PHI )+ZF*SIN( PHI ) 
XPORT = XPORTO+COS (THETA ) *COS(PHI ) +YPORTO*SIN ( THETA ) *COS (PHI > 

YOCAP = -XO*SIN(THETA)+YO*COS(THETA) 

YFCAP = -XF*SIN ( THETA )+YF*C0S( THETA) 

YPORT = -XPORTO*SIN(THETA)+YPORTO*COS(THETAJ i<ne 1 ajf 1 

ZOCAP = -XO»COS(THETA)^SIN(PHI)-YO*SIN(THETA)*SIN(PHI) +ZO*COS ( PHI ) 
ZFCAP = »X F*COS ( THETA) *SIN( PHI ) -YF*SIN( THETA) *SlN( PHI ) +ZF*C0S ( PHI) 

ZPORT = -XPORTO*COS ( THETA ) *SlN( PHI ) *YP0RT0*SIN( THETA ) *SIN (PHI ) 


C 

c 


A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 


10 

20 

30 

40 

50 

60 

70 

00 

90 

100 

110 

120 

130 

140 

150 

160 

170 

100 

190 

200 

210 

220 

230 

240 

250 

260 

270 

260 

290 

300 

310 

320 

330 

340 

350 

360 

370 

300 

390 

400 

410 

420 

430 

440 

450 

460 

470 

400 

490 

500 

510 

520 

530 

540 

550 

560 

570 
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60 


65 


70 


75 


60 


W 

to 


65 


90 


95 


100 


105 


110 


C . 

C . START OPTIMIZATION ♦ 

c . • 

c 

c 

INDEX s 0 

OLXCAP s ( XFCAP-XOCAPJ/50. 

C 

C 

c ♦ • 

C • FIRST FIND A CURVE WHICH FORCES THE HEADING OF THE • 

C • AIRCRAFT TOWARD T*E RUNWAY AT THE FINAL POINT • 

C . 

C . 

c 

slope = <yfcap-ypq:t)/(Xfcap-xport) 

YCURVE(l) e YOCAP 
AOT(i) = o. 

ADD Y ( 1 ) = 0. 

XCAP = XOCAP 
DO 10 I = 1 , 50 

XCAP s XCAP^DLXCAP 

EXPO s -5.*(XCAH-XFCAP)/(X0CAP-XFCAP) 

YCURVEU + 1) = (S| OPE* ( XCAP-XPORT) ♦( YPORT-YOCAP) UEXP( EXPO) +YOCAP 

AOY( 1*1) = -5./(X0CAp-XFCAP)*( YCURVE(I*1 J-YOCAp)* (SLOPE)* EXP (EXP 

1 0) 

ADD Y ( I +1 > = ( (-5,/{X0CAP-XFCAP) )**2)* ( YCURVE( 1*1 ) -YOCAP ) ♦ ( -5 ./ I X 
1 OCAP-XFCAP) >*SLOPE*EXP(EXPO)*2, 

10 CONTINUE 
COUNT = 0, 

c 

c 

C . 

C . INITIAL COST * 

C . 

c 

C 

XMXN = -40000 
XINC = 2500 
YMIN = -40000 
YINC = 2500 
H = p * 07 

CALL COST ( 0, 1. XOCAP. YOCAP •ZOCAP,XFCAP. OLXCAP, THETA t PHI* TOT AL.PNAL 
ITT) 

COST1 = TOTAL 
A = COSTl-PNALTY 

WRITE (6.9150) COUNT . COST1 . A t PNALTY 
WRITE (6.9220) 

WRITE (6.9230) ( (POSIT ( I , J) • J=1 . 3) ,1=1.51) 

WRITE (6,9030) 

00 20 I = 1.5 

WRITE (6.9040) I . ALFA (I) , BETA ( I ) 

20 CONTINUE 

WRITE (6,9050) 

30 DO 40 I = 1,5 

OALFA(I) = All 
OBETA(I) = A12 


A 

560 

A 

590 

A 

600 

A 

610 

A 

620 

A 

650 

A 

640 

A 

650 

A 

660 

A 

670 

A 

660 

A 

690 

A 

700 

A 

710 

A 

720 

A 

730 

A 

740 

A 

750 

A 

760 

A 

770 

A 

780 

A 

790 

A 

600 

A 

810 

A 

820 

A 

630 

A 

640 

A 

850 

A 

860 

A 

870 

A 

860 

A 

890 

A 

900 

A 

910 

A 

920 

A 

930 

A 

940 

A 

950 

A 

960 

A 

970 

A 

980 

A 

990 

A 

1000 

A 

1010 

A 

1020 

A 

1030 

A 

1040 

A 

1050 

A 

1060 

A 

1070 

A 

10B0 

A 

1090 

A 

1100 

A 

1110 

A 

1120 

A 

1130 

A 

1140 
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115 


120 


125 


130 


135 


1 H 0 

w 

I 

l jO 

1H5 


150 


155 


160 


165 


170 


C 

c 

c . • 

c . calculate gradient • 

c . • 

c 

c 

50 00 60 I = 1*5 

ALFA ( I ) = ALFACl)iDALFA(I) 

CALL COST ( 1 .0 , XOCAP. YOC AP t ZOCAP . XFCAP * OLXCAP • THETA « PHI . TOT AL * PN 
1 ALTY) 

COST 2 = TOTAL 

G Y ( I ) = (COST2-COST1)/ABS(UALFA(I)) 

IF (INDEX. EQ.O) AGY(I) = GT(I) 

IF (INDEX. EQ.l) BGY(I) = G Y ( I > 

WRITE (6.9160) I . GY ( I ) 

60 ALFA(I) s ALFACI)-DALFA(I) 

DO 70 I = 1.5 

BETA(I) = BETA ( 1 1 +DBET A ( I ) 

CALL COST (1 .O.XOCAP, YOCAP.ZOCAP.XFC AP.DLXCAP. THETA. PHI.TOTALiHN 
1 ALTY) 

C0ST2 = TOTAL 

GZ(I) = (COST2-CoSTl)/ABS(OBETA(I)) 

GZ(I) = 0. 

IF (INDEX. EQ.O) AGZ(I) = GZ(I) 

IF (INDEX. EQ.l) BGZ(I) = GZ 1 1 ) 

WRITE (6.9170) I.GZtl) 

70 BET A ( I ) = BETA ( I ) »DBETA < I ) 

IF (INDEX, EQ.l) GO TO 190 
GYMAX c ABS (GY (1 ) ) 

GZ1AX = ABS(GZd) ) 

DO 00 I = 2.5 

IF (GYMAX. LT, ABS(GY(I )) ) GYMAX s ABS ( G Y ( I ) ) 

80 IF (GZMAX.LT.ABS(G2(I))1 GZMAX = ABS(GZ(I)I 
C 

C 

c . 

C . DETERMINE SIZE OF STEP CHANGE • 

C . 

C 

C 

YALLOW s ( Y ALLOW- All )*0«95+All 
ZALLOW s ( Z ALLOW- A12) *0.95*A12 
IF (GYMAX. EQ.O.) YrATIO = 0. 

IF (GYMAX. NE.O.) YRATIO = YALLOW/GYMAX 

IF (GZMAX.EQ.O.) ZrATIO = 0. 

IF (GZMAX.NE.O.) ZRATIO = ZALLOW/GZMAX 

DO 90 I = 1.5 

ALFAOD(I) = ALFA ( I ) 

ALFA(I) c ALFA(I)-YRATIO*GY(I) 

BETAOO(I) = BETA ( I ) 

90 BET A ( I ) = BETA(I)-ZRATIO*GZ(I) 

CALL COST ( 0.0 .XOCAP.YOCAP. ZOCAP.XFC AP.DLXCAP, THETA «PHI .TOTAL. PNAL 
IT Y ) 

C0ST2 = TOTAL 

IF ( C0ST2 , GE , COST1 ) GO TO 150 
100 PRC ENT = ABS(C0ST2-C0ST1 )/COSTl 


A 1150 
A 1160 
A 1170 
A 1160 
A 1190 
A 1200 
A 1210 
A 1220 
A 1230 
A 12H0 
A 1250 
A 1260 
A 1270 
A 1260 
A 1290 
A 1300 
A 1310 
A 1320 
A 1330 
A 13H0 
A 1350 
A 1360 
A 1370 
A 1360 
A 1390 
A 1HO 0 
A 1 H 1 0 
A 1H20 
A 1H 30 
A 1HH0 
A 1H50 
A 1 H 6 0 
A 1H70 
A 1460 
A 1490 
A 1500 
A 1510 
A 1520 
A 1530 
A 1540 
A 1550 
A 1560 
A 1570 
A 1580 
A 1590 
A 1600 
A 1610 
A 1620 
A 1630 
A 1640 
A 1650 
A 1660 
A 1670 
A 1660 
A 1690 
A 1700 
A 1710 
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175 


180 


185 


190 


195 


200 


205 


210 


215 


220 


225 


C A 1720 

C A 1730 

C * . A 1740 

C . STOP CRITERION -- PERCENTAGE CHANGE IN COST INSIGNIFICANT * A 1750 

C . .A 1760 

C A 1770 

C A 1780 

IF (PRCCNT.GE.l.E-5) GO TO HO A 1790 

COUNT s COUNT + 1 A 1800 

CALL COST (Oil*XOCAP,YOCAPi20CAP*XFCAP*DLXCAP,THETA*PHI*TOTALiPNAL A 1810 

ITT) A 1820 

WRITE (6*9180) COUNT A 1830 

CALL MONIT ( COUNT t C0ST2 * PNALTY ) A 1840 

STOP A 1850 

110 CALL COST (0.1.XOCAP,YOCAP*ZOCAP*XFCAP*OLXCAP*THETAiPHIiTOTAL*PNAL A 1860 

1 TY I A 1870 

C0ST1 s TOTAL A I860 

COUNT = COUNT+1 A 1890 

A = C0ST1-PNALTY A 1900 

WRITE (6.9150) COUnT.COSTI * A.PNALTY A 1910 

WRITE (6.9220) A 1920 

WRITE (6.9230) ( ( POSIT ( I , J) , J=1 . 3 ) . 1=1 , 51 ) A 1930 

DO 120 I s 1.5 A 1940 

C A 1950 

C A I960 

C . • A 1970 

C . STOP CRITERION -- ALL GRADIENT COMPONENTS EQUAL TO ZERO . A I960 

C . .A 1990 

C A 2000 

C A 2010 

IF ( GY ( I ) • NE • 0 . ) GO TO 130 A 2020 

IF (GZ(I).NE.O.) GO TO 130 A 2030 

120 CONTINUE A 2040 

WRITE (6.9060) COUNT A 2030 

CALL MONIT (COUNT.COSTI. PNALTY) A 2060 

STOP A 2070 

130 WRITE (6.9070) A 2080 

DO 140 I = 1,5 A 2090 

WRITE (6,9190) I , ALF A ( I ) , BETA ( I ) A 2100 

140 CONTINUE A 2110 

C0ST2 = TOTAL A 2120 

C A 2130 

C A 2140 

C . .A 2150 

C , STOP CRITERION -- MAXIMUM NUMBER OF ITERATIONS REACHED . A 2160 

C . .A 2170 

C A 2180 

C A 2190 

IF (COUNT. LT.MAXIT) GO TO 30 A 2200 

WRITE (6,9200) A 2210 

CALL MONIT ( COUNT . COSTl , PNALTY ) A 2220 

STOP A 2230 

150 HALF =1 A 2240 

C A 2250 

C A 2260 

C . , A 2270 

C , REOUCE SIZE OF STEP CHANGE BT HALF • A 2280 
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230 

C 

C 

c 

• 

• 

IP COST HAS NOT DECREASED 

• A 

* A 

2290 

2300 

2310 


c 



A 

2320 




00 170 J = 1*3 

A 

2330 




00 160 I = 1 « 5 

A 

2340 

235 



ALFA ( I ) e (ALFA(I)+ALFAOD(I) )/2, 

A 

2350 



160 

BETA ( I ) = (RETA(I)+BETAOD(I) )/?. 

A 

2360 




half = J 

A 

2370 




WRITE (6*9210) HALF 

A 

2380 




CALL COST (0.0*XoCAP,YOCAP*ZOCAP,XFCAP.OLXCAP*THETAtPHI*TOTAL*PN A 

2390 

240 


1 ALTY) 

A 

2400 




C0SI2 = TOTAL 

A 

2410 




IF (C0ST2.LT. C0ST1) GO TO 100 

A 

2420 



170 

CONTINUE 

A 

2430 




HALF = 4 

A 

2440 

245 



inuex = 1 

A 

2450 




00 180 I s 1 i 5 

A 

2460 




DALFA(I) c -DALFAU) 

A 

2470 



180 

OBETA(I) s -DBETA(I) 

A 

2480 


c 



A 

2490 

250 

c 




2500 


c 

• 


• A 

2510 


c 

• 

PERTURB CURVE IN ThE OPPOSITE DIRECTION 

• A 

2320 


c 

• 


• A 

2530 


c 




2340 

255 

c 



A 

2550 




GO TO 50 

A 

2360 



190 

DO 200 I = 1*5 

A 

2370 




IF (AGY(I).LT.O.) GO TO 220 

A 

2380 




IF ( BGY ( I ) * LT * 0 « ) GO TO 220 

A 

2590 

260 



IF ( AGZ(I ) ,LT.O. ) GO TO 220 

A 

2600 




IF ( BGZ ( I ) . LT • 0 • ) GO TO 220 

A 

2610 



200 

CONTINUE 

A 

2620 




WRITE (6*9080) 

A 

2630 




DO 210 I c 1,5 

A 

2640 

265 



ALFA ( I ) e ALFAOD(I) 

A 

2630 



210 

BETA(I) e BETAOD(I) 

A 

2660 




CALL COST (0*1.XO^AP.YoCAP*ZOCAP*XFCAP*DLXCAP,THETA*PHIiTOTAL*PNAL A 

2670 



ITT) 

A 

2680 




CALL, MONIT < COUNT * COST 1 . PNALTY 1 

A 

2690 

270 



STOP 

A 

2700 



220 

8GYMAX c ABS (BGY ( 1 ) ) 

A 

2710 




BGZMAX c ABS(BGZd) ) 

A 

2/20 




DO 230 1 s 2*5 

A 

2/30 




IF (BGYMAX.LT.ABS(BGY(I))) BGYMAX r ABS ( BGY ( I ) ) 

A 

2/40 

275 


230 

IF (BGZMAX, LT.ABS(BGZd))) BGZMAX = ABS ( BGZ < I ) ) 

A 

2/50 



240 

WRITE (6*9210) HALF 

A 

2/60 


c 



A 

2770 


c 




2/80 


c 

• 


• A 

2790 

280 

c 

• 

CHECK EACH GRADIENT COMPONENT TO DETERMINE SIZE 

• A 

2800 


c 

• 

OF STEP CHANGE 

, A 

2810 


c 

• 


, A 

2820 


c 




2830 


c 



A 

2840 

285 



DO 320 I = 1*5 

A 

2830 
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IF ( HALF , EQ , 7 ) GO TO 250 

A 

2860 



IF (GYMAX.NE.O, ) AY = YALLOW/GYMAX/FLO AT ( HALF-3) *AGY ( I > 

A 

2870 



IF (BGYMAX.NE.O. ) BY s YALLOW/BGYM AX/FLOAT ( HALF-3 > *BGY t I ) 

A 

2860 



IF (GZMAX.EQ.O. ) AZ c 0. 

A 

2890 

290 


IF ( BGZMAX ,EQ* 0 • ) BZ * 0. 

A 

2900 



IF (GYMAX.EO.O. ) AY = 0. 

A 

2910 



IF ( BGYM AX « EQ * 0 • ) BY * 0. 

A 

2920 



IF (GZMAX,NE.O, ) AZ e Z ALLOW/GZMAX/FLOAT ( HALF-3 > *AGZ ( I ) 

A 

2930 



IF (BGZMAX. NE.O.) BZ e ZALLOW/BGZMAX/FLOAT ( HALF-3 > *BGZ ( I ) 

A 

2940 

295 


GO TO 260 

A 

2950 


250 

AY s -DALFA(I) 

A 

2960 



BY c DALFA(I) 

A 

2970 



AZ = -DBETA(I) 

A 

2980 



BZ = DBETA(I) 

A 

2990 

300 

260 

IF (AGY(I).LE.O.) GO TO 270 

A 

3000 



IF (BGY(I).GE.O.) ALFA(I) = ALFAOD(I) 

A 

3010 



IF (BGY(I).LT.O.) ALFA ( I ) = ALFAOD(I)+BY 

A 

3020 



GO TO 290 

A 

3030 


270 

IF (AGY(I).LT.O.) GO TO 280 

A 

3040 

305 


IF ( BGY ( I ) ,GE, 0* ) ALFA(I) * ALFAOO(I) 

A 

3050 



IF (BGY(I).LT.O.) ALFA(I) = ALFAOD ( I ) +BY 

A 

3U60 



GO TO 290 

A 

3070 


280 

IF (AGY(I)'LT.BGy(I) ) ALFA(I) s ALFAOO(I)-AY 

A 

3060 



IF (AGY(I).GE.BGY(I)) ALFA(I) s ALFAOD(I)+BY 

A 

3090 

310 

290 

IF (AGZ(I).LE.Oo GO TO 300 

A 

3100 



IF ( BGZ ( I ) « GE • 0 « ) BETA(I) = BETAOD(I) 

A 

3110 



IF (BGZ(I).LT.O.) BETA(I) = BETAOO(I)+BZ 

A 

3120 



GO TO 320 

A 

3130 


300 

IF < AGZU ) ,LT,0. ) GO TO 310 

A 

3140 

315 


IF ( BGZ { I ) » GE« 0 • ) BETA ( I ) « BET AOD ( I ) 

A 

3150 



IF (BGZ(I).LT.O.) BETA(I) * BETAOO(I)+BZ 

A 

3160 



GO TO 320 

A 

3170 


310 

BET A ( I ) = BET A ( I ) - AZ 

A 

3180 


320 

CONTINUE 

A 

3190 

320 


CALL COST ( 0. 0 t XOC AP, YOCAP. ZOCAP • XFC AP. DLXCAPt THETA* PHI ♦ TOT AL.PNAL 

A 

3200 


ITT) 

A 

3210 



COST2 = TOTAL 

A 

3220 



IF (C0ST2.LT.C0ST1) GO TO 420 

A 

3230 



HALF r HALF+1 

A 

3240 

325 


IF ( HALF , LT • 7 ) GO TO 240 

A 

3250 



WRITE (6*9210) HALF 

A 

3260 

1 


GYMIN s AGY(l) 

A 

3270 



J = 1 

A 

3280 



GZMIN c AGZ(l) 

A 

3290 

330 


K « 1 

A 

3300 



00 340 I = 2*5 

A 

3310 



IF (GYMIN, LE. AGY(I) ) GO TO 330 

A 

3320 



GYMIN b AGY(I) 

A 

3330 



J = I 

A 

3340 

335 

330 

IF (GZMIN, LE, AGZ(I) ) GO TO 340 

A 

3330 



GZMIN = AGZ(I) 

A 

3360 



K s I 

A 

3370 


340 

CONTINUE 

A 

3380 



00 360 I = 1*5 

A 

3390 

340 


IF (GYMIN. LE,BGY(IM GO TO 350 

A 

3400 



GYMIN b BGY ( I ) 

A 

3410 



J * 1 + 5 

A 

3420 
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W 

1 


345 


350 


355 


360 


365 


370 


375 


380 


385 


390 


395 


350 IF <GZMIN.LE*BGZ{I>) GO TO 360 
GZMIN = 8GZ ( I ) 

K = 1 + 5 

360 CONTINUE 

IF ((GYMIN.LT. 0.0). OR, (GZMIN.LT. 0*0)) GO TO 370 

CALL COST ( 0. 1 * XOCftP. YOCAP , ZOC AP , XFC AP , DLXCAP , THETA* PHI* TOT AL.PNAL 
1TY ) 

COUNT * COUNT+1 

WRITE (6*9090) COUNT 

CALL MONIT ( COUNT ,COSTl » PNALTY ) 

STOP 

370 DO 380 I = 1.5 

ALFA ( I ) = ALFAODfl) 

380 BETA(I) = BETAOD(I) 

IF ((GYMIN.LT,0.0) # ANO. (GZMIN. Gt. 0.0)) GO TO 390 
IF ( (GYMIN.LT.O. ) .AND. (GZMIN.LT. 0. )) GO TO 400 
IF (K.LC.5) BETA(K) s BETA(K)-DBETA(K) 

IF (K.GT.5) BETA(K-5) c BETA ( K-5 ) +OBETA ( K-5 ) 

GO TO 420 

390 IF (J.LE.5) ALFA(J) = ALF A ( J) -DALFA ( J ) 

IF (J.GT.5) ALFA ( J-5 ) = ALFA ( J-5 ) +DALF A ( J-5 ) 

GO TO 420 

400 IF (J.LC.5) ALFA ( J ) c ALFA ( J ) -DALF A ( J) 

IF (J.GT.5) ALFA ( J-5 ) c ALFA ( J-5 ) +DALFA ( J-5 ) 

IF (K.LC.5) BETA(K) = BCTA(K>-OBETA(K> 

IF (K.GT.5) BET A ( K-5 ) = BETA ( K-5 ) +DBET A ( K-5 ) 

CALL COST ( 0.0. XOCAP.YOCAP. ZOCAP.XFCAP.DLX CAP, THETA .PHI. TOTAL, PNAL 
1TY ) 

C0ST2 = TOTAL 

IF (C0ST2.LT.C0ST1) GO TO 420 
DO 410 I r 1,5 
410 BET A ( I ) = BETAOD,I) 

420 INDEX = 0 

CALL COST ( 0,1. XOCAP.YOCAP, ZOCAP.XFCAP, OLXCAP, THETA, PHI, TOTAL, PNAL 
1TY) 

GO TO 100 
C 
C 

9010 FORMAT ( 5X , 14HINI TaL X.Y.Z5 • 3 ( F12 , 2 , 3X ) , 7H MET ErS , / , 5X , 13HF I N AU X 
1.Y.ZJ , 3 ( F12 « 2, 3X ) « 7H METERS */. 5X , 23HAIRP0RT LOCATION, X.Yl ,2(F12 
2 . 2 * 3X ) , 7H METERS) 

9020 FORMAT ( 5x , 43HPERTURB TRAJECTORY IN Y AND Z DIRECTONS BY ,F6.2f5H 
1 ANO , F6 • 2 , 4 2H METERS, RESPECTIVELY FOR CALCULATING GRAD, 5HIENTS ) 
9030 FORMAT ( 1 3X , 4HALF A , 16X , 4HBET A ) 

9040 FORMAT ( 1 0 X , I 1 , 1PE16 .9 , 4X , 1PE16 . 9) 

9050 FORMAT (////) 

9060 FORMAT (////, IX , 13HAT ITERATION .I2.49H ALL GRADIENTS EQUAL TO 
1 ZERO, PROGRAM STOPS) 

9070 FORMAT ( 10 X , 2HNO , IX . 4HALF A , 16X , 4HBET A ) 

9080 FORMAT (5X.43HALL GRADIENTS PERTURBED BOTH DIRECTIONS > 0) 

9090 FORMAT (1X.13HAT ITERATION tl2,16H OPTIMUM REACHED) 

9100 FORMAT ( 3 A 1 0 , / , 4 A l n ) 

9110 FORMAT ( 1H1 » 20X , 3 AlO , // , 4A10 , //// ) 

9120 FORMAT ( IX . 19HINF0AMATION INPUTS . // , 5X , 21HMAXIMUM ITERATION SET.l 
1HS.I3) 

9130 FORMAT ( 5X , 47HMAXIMUM ALLOWED CHANGES PER ITERATON IN Y AND Z,27H 
10IRECTI0NS. RESPECTIVELY* .1PE10.3.5H ANO ,1PE10.3,7H METERS) 


A 3430 
A 3440 
A 3450 
A 3460 
A 3470 
A 3480 
A 3490 
A 3500 
A 3510 
A 3520 
A 3530 
A 3540 
A 3550 
A 3560 
A 3570 
A 3580 
A 3590 
A 3600 
A 3610 
A 3620 
A 3630 
A 3640 
A 3650 
A 3660 
A 3670 
A 3680 
A 3690 
A 3700 
A 3710 
A 3720 
A 3730 
A 3740 
A 3750 
A 3760 
A 3770 
A 3780 
A 3790 
A 3800 
A 3810 
A 3820 
A 3830 
A 3840 
A 3850 
A 3060 
A 3870 
A 3880 
A 3890 
A 3900 
A 3910 
A 3920 
A 3930 
A 3940 
A 3950 
A 3960 
A 3970 
A 3980 
A 3990 
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*+00 9140 FORMAT <5X#22HINITIAL ALFA AND BETA ! , / 1 13X < 4HALFA * lbX # 4HBET A # 5 ( / f 1 A 4000 

10X»U#1X« 1PE16.9,4X#1PE16.9M A 4010 

9150 FORMAT t //// 1 IX ♦ IOhITER AT I ON . 1 3 , // , 5X t 1 4HT0T AL COST IS ,1PE16.9#/ A H020 
1 ♦ 5X t 22HTRUE ANNOYACE(NII) IS i IPE 16 . 9 # / * 5X * 42HPEN ALT Y DUE TO Al A 4030 

2RCRAFT CONSTRAINTS IS .1PE16.9#//) A 4040 

^05 9160 FORMAT ( 1 OX # 12 * 17HrH Y-GRADIENT IS #1PE16,9) A 4050 

9170 FORMAT < 10X , 12 , 17HTH Z»GR ADIENT IS .1PE16.9) A 4060 

9160 FORMAT <////tlX.13HAT ITERATION #I2#24H PERCENTAGE CHANGE IN CQ*33 A 4070 

1 HST LESS THAN .001*# PROGRAM STOPS) A 4000 

9190 FORMAT ( 1 OX # I 1 # 2X # 1PE16 # 9 # 4X # 1PE1 6 . 9 ) A 4090 

^10 9200 FORMAT ( 10 X # 41HREACH MAXIMUM ITERATION SET. PROGRAM STOP) A 4100 

9210 FORMAT ( 10X # 7HHALF = ,12) A 4110 

9220 FORMAT ( 1 OX , 10HTRA JECTORY , / ♦ 10X # 12HX COORDINATE » 0 X , 12HY C00RDINA1E A 4120 
1 # flX # 12HZ COORDINATE#/. 12X • 7H C METER )# 13X # 7H< METER )# l3X . 7H ( METER M A 4130 
9230 FORMAT (10X.3 C 1PE16.9.4X) ) A 4140 

415 END A 4150 


45000B CM STORAGE USED 7.620 SECONDS 


bd 

l 

oo 
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SUBROUTINE 


5 


10 


15 


20 


25 


30 


35 


HO 


45 


50 


55 


I 

i 


COST 73/172 TS 


FTN 4.6+452 04/27/79 11,45,47 PAGE 1 


SUBROUTINE COST ( IGRAO, I WRITE. XOCAP, YOCAP, ZOCAP.XFCAP.OLXCAP, THETA 
liPHI .TOTALtPNAUTY) 

COMMON ALFA ( 5 ) . BET A ( 5 ) » POSI T < 53 , 3 ) , ARR AY t 578 . 9 ) % NMAP 
COMMON /CURVE/ YCURVE ( 51 ) » ADY ( 51 ) , AOOY ( 51 ) 

COMMON /AIRPORT/ XPORT , YPORT * ZPORT 
COMMON /AC/ XtY.Z 
EXTERNAL FCN 
PNALTY = 0. 

XCAP = XOCAP 
PI = ATAN ( 1 , ) *4. 

C2 = Pl/ABS(XFCAP-XOCAP) 

C3 s ABS<XFCAP-XOCAP)/4. 

DO 10 I = ItNMAP 
ARRAY (1,4) s 0. 

10 ARRAY (1*5) s 0. 

C 

C *•••• 

C . 

C » MULTIPLY BY EXPONENTIAL TERM SUCH THAT THE FINAL • 

C . HEADING OF AIRCRAFT IS TOWARD THE RUNWAY 
C , 

C 

C 

00 50 I = 1.51 

Y 2 = l.-EXP(-(XFcAP.XCAP)/C3) 

Y 5 s (Y2-l,)/C3 

Y 9 s 0.0 

Y 8 = Y9 
Y7 s Y6 
Y 6 r Y7 
Y3 = Y6 

C 

c 

C . • 

c . GENERATE SINE HARMnNICS • 

C . 

C • 

c 

DO 20 J = 1*5 

TRIGOX = FLOAT ( J )♦( XCAP -XOCAP ) *C2 
Y3 s Y3+ALFA< J)*SIN(TRIGOX) 

Yfl s YQ+BETAt J) ♦SlNURIGOX ) 

Y6 5 Y6 + FLOATU>*C2*ALFA< J) *COS( TRIGOX ) 

Y7 = Y7-FLOAT( J**2)*(C2**2)*ALFA( J)*S1N(TRIG0X) 

20 Y9 = Y9+FLOAT(J)*C2+BETA(J)*COS<TRIGOX) 

OLYCAP s Y2*Y3 
0L2CAP r Y2*YB 
ZCAP = zocap+dlzcap 
YCAP = DLYCAP+YCURVE(I) 

C 

c 

C . 

c . AIRCRAFT CONSTRAINTS • 

C . 

C 

C 


B 10 
B 20 
B 30 
B 40 
6 50 

B 60 
B 70 
B 80 
B 90 
B 100 
B 110 
B 120 
B 130 
B 140 
B 150 
B 160 
B 170 
6 180 
B 190 
B 200 
B 210 
B 220 
B 230 
B 240 
B 250 
B 260 
B 270 
B 280 
B 290 
B 300 
B 310 
B 320 
B 330 
B 340 
B 350 
B 360 
B 370 
B 380 
B 390 
B 400 
B 410 
B 420 
B 430 
B 440 
B 450 
B 460 
B 470 
B 400 
B 490 
B 500 
B 510 
B 520 
B 530 
B 540 
B 550 
B 560 
B 570 


OY » Y2*Y6+Y3*Y5 
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60 


65 


70 


75 


60 


65 


90 


95 


f 

100 


105 


110 


DY = OY+ADY ( I ) 

ODY = Y2*Y7*2.*Y5*Y6+Y3*Y5/C3 
DOY = DDY-f AOOY ( I ) 

DDY = OOY/(1+DT**2) 

OZ s Y2*Y9+Y5*Y8 
DZ = OZ+TANtPHI) 

OZ = 0. 

PNALTY = PNAUTY+(DDY/.001)**(20)+(DZ/.14)**<20) 

X = XCAP*COS(THETA)*COS(PHl )-YCAP*SI Nt THETA >-ZCAP*COS( THETA )*SIN 

1 (PHI) 

Y = XCAp*SIN(THETA)*COS(PHl )+YCAP*C0S (THETA )«ZCAp+SIN( THETA )*SIN 
1 (PHI) 

Z = XCAP*SIN(PHI)+ZCAP*COS(PHI) 

DO 40 K s ItNHAP 

RANGE s ((X-ARRAY(Ktl)>**2+(Y-ARRAY(K,2))**2+Z**2)**.5 
DB = 115.-22.5*ALOG10(3.2ai*RANGE/500.) 

IF (0B.LE.ARRAY(K*4) ) GO TO 4o 
ARRAY ( K . 4 ) = R 

IF ( ARRAY ( K i 4 ) .LT. 55# ) GO TO 40 
IF (ARRAY(K«3).EQ.0. ) GO TO 30 
C 

C 

c ♦ 

C • ANNOYANCE INTEGRATION OVER A SINGE BLOCK • 

C . 

C 

c 

SMALLP s ARRAY(K,3)/( ARR AY < K « 7 ) - ARR AY < K i 6 ) )/( ARRAY ( K . 9 ) -ARRAY ( 
1 K « 8 ) ) 

CALL GAUSS (ARRAY(K*6) * ARRAY(Kf 7) , ARRAY(KtB) »ARRAY(K.9 ) »FCNi TE 

1 MP) 

ARRAY ( K ♦ 5 ) c TEMP*SMALLP 
GO TO 40 

30 ARRAY ( K • 5 ) = 0. 

40 CONTINJE 

IF ( I WRITE .EQ. 0 ) GO TO 50 
II = I 

P0SIT(II,1) = X 
POSIT ( 1 1 « 2 ) = Y 
POSIT (11*3) = Z 
50 XCAP = XCAP+DLXCAP 
C 

C ♦••••••• 

C . 

C . TOTAL POPULATON EXPOSED TO NOISE ABOVE 55 EPNOB 
C . 

C 

c 

PEOPLE = o. 

DO 60 K s l.NMAP 

IF (ARRAY (Ki5 ) .Eq . 0 • 0 ) GO TO 60 
PEOPLE s ARRAY(K.3)+PEOPLE 
60 CONTINUE 
FX = 0. 

DO 70 K t 1 « NMAP 

ARRAY ( K » 5 ) = ARR AY ( K 1 5 ) /PEOPLE 
FX a FX + ARRAY ( K « 5 ) 


B 560 
B 590 
B 600 
B 610 
B 620 
B 650 
B 640 
B 650 
B 660 
B 670 
B 660 
B 690 
B 700 
B 710 
B 720 
B 730 
B 740 
B 750 
B 760 
B 770 
B 780 
B 790 
B 600 
B 810 
B 820 
B 630 
B 840 
B 850 
B 660 
B 870 
B 660 
B 890 
B 900 
B 910 
B 920 
6 930 

8 940 

B 950 
B 960 
B 970 
B 980 
B 990 
B 1000 
B 1010 
B 1020 
B 1030 
B 1040 
B 1050 
B 1060 
B 1070 
B 1080 
B 1090 
B 1100 
B 1110 
B 1120 
B 1130 
B 1140 
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115 70 CONTINUE 

total = fx*pnalty 

RETURN 

END 

4 1 00 OB CM STORAGE USEO .874 SECONOS 


8 1150 
B 1160 
B 1170 
B 1160 
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PAGE 1 


5 


10 


15 


20 


25 


30 


35 


40 


45 


subroutine: monit <ia,aa*bb) 

COMMON ALFA (5) t BET A (5) , POSIT (53*3) * ARR A Y ( 578 , 9 ) . NMAP 
COMMON /LABEL/ LINpO { 4 ) *LLOC ( 3 ) 

COMMON /SCALE/ XMIN , XINC • YMIN , YINC 
OIMENSION PCRIT(IO) 

DIMENSION XM(1026), YM(1026) 

DIMENSION XP ( 53 ) * YP(53>* 2P(53). NA(5)* NB(3) 

EQUIVALENCE (XMtl) •ARRAY(lil) ) * ( YM ( 1 ) * ARRAY < 1 , 2 I ) 

EQUIVALENCE (XP(1) *POSIT(l*l) ) * ( YP ( 1 > * POSI T< 1 • 2 ) ) * ( 2P ( 1 ) * POSIT ( 1 
1.3) ) 

DATA NB/lOHTOTAL POPU* 10HLATION ANN*9H0 YAnCE = / 

C 

C 

c . • 

c , DOCUMENTATION • 

C , 

c 

c 

cc = aa-bb 

WRITE (6*9010) 

WRITE (6*9020) XA«AA*CC«BB 
DO 10 X = 1*51 

WRITE (6*9030) (POSIT(I*J)*J=l*3> 

10 CONTINUE 

WRITE (6*9040) 

DO 20 I = 1 * NMAP 

WRITE (6*9050) ( ARRAY ( I * J ) « J*l • 5 ) 

20 CONTINUE 

WRITE (97*9060) ( ( POSIT ( I • J ) • J = 1 * 3 ) • 1 = 1 * 51 ) 

WRITE (97,9070) ( ( ARRAY ( I * J ) « J=1 • 5 ) *1 = 1* NMAP) 

RETURN 

C 

9010 FORMAT ( 10X *55H0PTIMUM TRAJECTORY FOR LANDING AT PATRICK HENRY AIR 
1PORT,//,10X*59HNOISE BELOW 55 EPNDB IS CONSIDERED NUT NOISY, ANNUy 
2ANCE = 0*/* 10X*23HUNIT FOR NOISE IS EPN03 * / * 1 OX * 50 HUNIT FOR COORUI 
3NATES* AIRCRAFT TRAJECTORY IS METER*//) 

9020 FORMAT (10X*l3HAT ITERATION • 12. /. 15X * 14HT0TAL COST IS *1PE16.9*/* 
1 15X , 23HTRUE ANNOYANCE ( Nl I ) IS « 1PE1 6, 9 * / * 15X , 1 1 HpENALTY IS *1PE16, 
29*//,lOX,lBHOPTlMUM TRAJECTORY*/* 10X, 12HX COORDINATE , 8X , 12HY COOKD 
3 1 NATE * 8X * 12H2 COORD I NATE •/* 12X* 7H ( METER )* 13X , 7H ( METER )* 13X • 7H ( METE 
4R) ) 

9030 FORMAT ( 10X * 3 ( 1PE16, 9. 4X ) ) 

9040 FORMAT ( // . 10 X * 32HPOPUL ATION-NOISE-ANNO YANCE CHART ,/• 10X • 10HX-POSI 
1TION,5X.10HY-POSITION*5X,10HPOP, INDEX. 5X. HHNOISE LEVEL • 4X • 9HANN0 
2YANCE) 

9050 FORMAT ( 10X , 3 (F10 , 3* 5X ) « 2 ( 1PE10 , 3. 5X) ) 

9060 FORMAT (3E12.6) 

9070 FORMAT (5E12.6) 

ENO 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


10 

20 

30 

40 

30 

60 

70 

80 

90 

100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

2S0 

290 

300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

4Q0 

490 


4100 OB CM STORAGE USED 


,284 SECONDS 
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5 


10 


15 


20 



SUBROUTINE GAUSS ( XN • XX • YN » YX , FCN , F INT ) 


U 

10 


COMMON /AC/ XAiYA»ZA 


0 

20 


DIMENSION X ( 5 ) • F(5)t XI ( 5 ) • W ( 5 ) 


0 

30 


DATA XI *W,N/-0. 577350269 i 0 , 577350 269 t 0 ♦ 0 i 0*1.*1* »0*Qt 0*2/ 


0 

40 




0 

50 




D 

60 



• 

0 

70 


GAUSSIAN QUADRATURE INTEGRATION WITH FOUR POINTS 

• 

D 

60 



* 

D 

90 




0 

100 




0 

110 


DO 10 I = 1 « N 


U 

120 


r ( I > * <YX-YN)/2.*XI(I)+(YX4YN)/2. 


U 

130 

10 

x<l> = (XX-XN)/2 # *XI(I)*(XX*XN)/2. 


0 

140 


FINT = 0. 


D 

150 


DO 30 J = liN 


D 

160 


F(J) = 0« 


0 

170 


DO 20 I s 1 • N 


0 

160 

20 

F ( U ) * F(J)+WtI)*FCN(X(I),Y(J)) 


U 

190 


F(J) = F ( J) * (XX**XN ) /2 • 


0 

200 

30 

FINT = FINT+W ( J) *F ( J ) 


D 

210 


FINT = FINT* ( YX-YN \ /? 


U 

220 


RETURN 


0 

230 


END 


0 

240 


PAGE 1 


410008 CM STORAGE USED 


166 SECONDS 
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FUNCTION FCN 73/172 TS 


FTN H#6+H52 OH/27/79 11.H5.H7 


5 


FUNCTION FCN <XtY> E 

COMMON /AC/ XA « YA < Z A £ 

RANGE! = SORT( (X-XA)**2*(Y-YAU*2 + ZA**2) l 

ARG s 129 , 12*22 *5* ALOGlO {RANGE ) t 

FCN = (3*36E«*6*10.**t ,103*ARG) )/( # 2*10«**< . 03* ARG > *1 . **3E-H* 1 O* ** < • £■ 

1 08*ARG ) ) E 

RETURN L 

END E 


10 

20 

30 

HO 

50 

60 

70 

80 


H1000B CM STORAGE USED 


•100 SECONDS 


PAGE 1 
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>» COST REPORT FOR LISTOAF <<< 

04/27/79 11,45.59 

RESOURCE BILLING RATf UNITS USED LOST 


CENTRAL PROCESSOR 

*105.00 

/HOUR 

9,314 

CP SECONDS $ 

.27 

PERIPHERAL PROCESSOR 

20.00 

/HOUR 

9.737 

PP seconds 

• 05 

I/O 

80.00 

/HOUR 

2*926 

io seconds 

.07 

FIELD LENGTH 

3,00 

/KILO-WRD-HOUR 

205.576 

KILO-JRO-SECS, 

.17 

(BASIC COST EXCLUDES LINES PRINTED, 
AND PLOTTER TIME CHARGES) 

CARDS PUNCHED 


basic COST 

.56 

JOB PRIORITY 3 

PRIORITY COST FACTOR 1.00 

APPROXIMATE AOJUSTEO COST 

.56 


AS OF LAST ACCOUNT UPDATE i ACCOUNT EXPIRES 04/30/79, FUNDS LEFT S 6037.31 


04/27/79 UVA NOS/BE 1,2 LEVEL 454-03/11/75 

11. 45. 47. LISTOAF FROM ♦GD/AB 

11.45.47. LIST,M3117A«T100. 

11.45.47. ATT ACHtQ i NEWTIOY. 

11. 45. 47. PF CYCLE NO. = 002 

11.45.47. FTN(I=Q) 

11.45.59. 450008 CM STORAGE USED 

11.45.59. 9.292 CP SECONDS COMPILATION TIME 

11.45.59. STOP 

11. 46, 00. EJ END OF JOB « A8 


PRINT COST $000.88 LISTOAF //// CND OF LIST //// 0000805 LINES 



DISTRIBUTION LIST 


Copy No . 
1-2 


3-5 


6-7 

8-9 

10 

11 

12 

13 - 15 


16 


NASA Scientific and Technical Information Facility 
P.0. Box 8757 

Baltimore/Washington International Airport 
Baltimore, HD 21240 

Mr. Richard DeLoach 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Va 23665 

I. D. Jacobson 

G . E . Cook 

L. S. Fletcher 

E. A. Parrish, Jr. 

I. A. Fischer 

Office of Sponsored Programs 
E. H. Pancake 

Science/Technology Information Center 
Clark Hall 

RLES Files 


9967:hfc 



UNIVERSITY OF VIRGINIA 


School of Engineering and Applied Science 

The University of Virginia's School of Engineering and Applied Science has an undergraduate 
enrollment of approximately 1,000 students with a graduate enrollment of 350 There are approximately 
120 faculty members, a majority of whom conduct research in addition to teaching 
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the doctorate, the Biomedical and Materials Science Departments grant only graduate degrees 

The School of Engineering and Applied Science is an integral part of the University (approximately 
1,400 full-time faculty with a total enrollment of about 14,000 full-time students), which also has 
professional schools of Architecture, Law, Medicine, Commerce, and Business Administration In addition, 
the College of Arts and Sciences houses departments of Mathematics, Physics, Chemistry and others 
relevant to the engineering research program This University community provides opportunities for 
interdisciplinary work in pursuit of the basic goals of education, research, and public service 
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